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Abstract 

We have applied the Mehiikov criterion to examine a global homoclinic bifurcation and 
transition to chaos in a case of the Duffing system with nonlinear fractional damping and 
external excitation. Using perturbation methods we have found a critical forcing amplitude 
above which the system may behave chaotically. The results have been verified by numeri- 
cal simulations using standard nonlinear tools as Poincare maps and a Lyapunov exponent. 
Above the critical Melnikov amplitude /ic, which is the sufficient condition of a global ho- 
moclinic bifurcation, we have observed the region with a transient chaotic motion. 
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1 Introduction 



A nonlinear oscillator with the Duffing term but linear damping is one of the simplest systems 
leading to chaotic motion studied originally by Ueda ^ [21 01 • The problem of its nonlinear 
vibrations has attracted researchers from various fields of research; from natural science [H [31 
and mathematics ()^ to mechanical El El and electrical engineering 0l2lEl- This system, 
for a negative linear part of stiffness, shows homoclinic orbits, and the transition to chaotic 
vibration can be treated analytically by the Melnikov method. Note that such a treatment 
has been already performed successfully to selected problems with various potentials [HI 13. In 
spite of fact that a large bibliography is devoted to vibration of a single realization [1-12] or 
set of coupled Duffing oscillators |E1 El El E| with numerous modifications to potential 
and forcing parts, the problem of nonlinear damping in chaotically vibrating system has been 
focused shortly [12j . The fractionally damped Duffing's equation, for a single well potential 
with positive linear part, has been studied in detail by Padovan and Sawicki |18j . They were 
especially focused on the determination of the influence of fractional damping on the frequency 
amplitude response. Some additional insight into this problem can be also found in the context 
of self excitation effects |T2 El El El E]- On the other hand, in the paper of Trueba et 
al. ;f2j . the systematic discussion on square and cubic damping effects on global homoclinic 
bifurcations in the Duffing system was given given while Awrejcewicz and Holicke applied the 
Melnikov theory to Duffing system in a dry friction limit j2Uj . In the present paper we reexamine 
the problem of nonlinear damping including fractional damping cases by both; the analytical 
Melnikov method as well as numerical simulations. Namely we are going to study a single degree 
of freedom nonlinear oscillator with the Duffing potential and fractional damping: 

X + ax \x\^~'^ + 5x + '^x^ = ^ cos ujt, (1) 

where x is displacement and x velocity, respectively, while the external potential force Fx and 
corresponding potential V{x) (Fig. ^) are defined; 

Here the non-linear damping term is defined by the exponent p: 

dpt(x) = Qx . (3) 

We have plotted the above function versus velocity {v = x) in Fig. ^ for few values of p. Note 
that, the case p — > (see p = 0.1 in Fig. ^ for a relatively small velocity) mimics the dry 
friction phenomenon [211 122j . 



2 Melnikov Analysis and Numerical Simulations 

We start our analysis from the unperturbed Hamiltonian H^: 

H'='^ + V{x). (4) 

The potential function V{x) (Fig. ^) has the local peak at the saddle point x = 0. Existence 
of this point with a horizontal tangent makes possible homoclinic bifurcations to take place. 
This includes transitions from regular to chaotic solutions. To study effects of damping and 
excitation, we apply small perturbations around the homoclinic orbits. This implies using a 



small parameter e in the Eq. ^ with perturbation terms. Equation ^ can be simultaneously 
uncoupled into two differential equations of the first order: 



X = V (5) 
V = — eav\v\^~^ — 6x — jx^ + ejl cos Lot, 

where ea = a and efl = fi, respectively. 

At the saddle point a; = 0, for an unperturbed system (Fig. the system velocity reaches 
zero in velocity v = (for infinite time t = ±00) so the total energy has only its potential part. 
Transforming Eqs. 2,4 for a chosen nodal energy (E = 0) and for (5 < 0, 7 > we get the 
following expression for velocity: 



dx /„ / fe^ 7X' 



dt y V 2 

Now one can perform integration over x: 



4" 

(6) 



t-to = ± I , (7) 

2 

where to represents an integration constant. Finally, we get so called homoclinic orbits (Fig. |2j: 
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X 



7 cosh (\/^(t - to)) 

7 cosh {^{t - to)) 

where '+' and '— ' signs are related to left- and right-side orbits, respectively (Fig. |2j. Note, 
the central saddle point xo = is reached in time t corresponding to +00 and —00, respectively. 

In case of perturbed orbits (a stable manifold) and (an unstable manifold) the 
distance between them is given by the Melnikov function M(to): 

/+00 
h{x*,v*) Ag{x*,v*)dt (9) 
-co 

where the corresponding differential form h means the gradient of unperturbed Hamiltonian 
(Eq. EI): 

h = {6x* + 7(x*)3) dx + v*dv, (10) 
while 5 is a perturbation form (Eq. ^ to the same Hamiltonian: 

g = (fl cos ujt — av* \v*f^^] dx. (11) 



All differential forms are defined on homoclinic orbits {x,v) = {x*,v*) (Eq. [H)). 

Condition for a global homoclinic transition, corresponding to a horse-shoe type of stable and 
unstable manifolds cross-section (Fig. can be written as: 

V M(to) = and / 0. (12) 

Evaluating the corresponding integral (Eq. ^ after some lengthly algebra the last condition 
(Eq. I12|l yields to a critical value of excitation amplitude |^^] Hc- 

JZ!lt^B(l±^,l±l)c..J^], (13) 



where the B(r, s) is the Euler Beta function: 



B(r,s) 



r(r)r(g) 

r(r + s) 



(14) 



and r(n) denotes the Euler Gamma function: 



r(n + 1) = nr(n). 



(15) 



In Fig. I2K we plotted the results of Melnikov analysis for a critical amplitude fj-c/o: for few 
values of p while 6 = —1, ■j = 1. For fj, > Hc the system can transit to chaotic vibrations. In 
the next figure (Fig. IHh) we plotted /ic versus p. Here we have also compared the results of 
analytical investigations represented by a solid line with numerical simulations represented by 
circle points. The numerical results can give a local criterion of chaotic vibration appearance, 
namely a positive value of the largest Lyapunov exponent Ai Note that in all cases (Fig. 
ISb) points are above the line. This can be associate with the fact that the Lyapunov exponent is 
a sufficient criterion of system transition to chaos while the Melnikov condition is only necessary 
one. Thus we can conclude that both sets of results: analytical and numerical are in good 
agreement. The examples of the largest Lyapunov exponent Ai versus the excitation amplitude 
fj, are plotted in Figs. \^ and b. In our numerical code we started calculations from the same 
initial conditions (xo,wo) = (0.45,0.1) for every new value of /i. One can clearly see the point 
of Ai sign change. This indicates transition to chaotic vibrations. Note that for p = 0.1 we 
have got some small but positive values of Ai for weak excitation (/i close to 0). When we 
examined this case with more details we observed that amplitude of vibrations was very small. 
Appearance of small positive value of Ai seams to be a combined effect of physical phenomenon 
of dry friction and our assumed numerical integration procedure of constant time steps At (in 
most cases we used At = 27r/(500u;)) . It is known from the studies of a stick-slip phenomenon 
j24| I25j that in such cases integration should be performed very carefully. Indeed, making time 
steps relatively smaller (50 times smaller) we checked that it is sufficient to reduce the positive 
value around /i = leaving, in principle, unchanged that Ai which appears at ~ 0.27. To 
reduce this artificial effect of small positive values of Ai and assuming an efficient algorithm (not 
so much time consuming) we decided to put a minimal threshold of Ai = 0.1 above which the 
system was identified in a chaotic state (see Figs. ISbiEt^ and^). 

For a given value p (p = 0.1) we show (in Fig. for ^ = 0.27 and Fig. Eh for /i = 0.25) 
typical phase portraits (represented by lines) and simultaneously Poincare maps (represented by 
points). Note that in Fig. EH we plotted trajectories for t G [1200, 1500], in arbitrary units of l/cu, 
while for Fig. Eb we found transient behaviour. To show differences we plotted corresponding 
trajectories for t G [0, 1200] and t G [1200, 1500] in various colours. The characteristic time 
t = 1200 has been chosen in such a way to reach by our system a steady state. It is easy 
to see that Fig. [5^ represents a non-periodic vibration state (/i = 0.27) just after the system 
transition to a chaotic motion (see function Ai versus /x in Fig. 0^). The corresponding Lyapunow 
exponent Ai = 0.175 is positive indicating on exponential divergence. In spite of a small number 
of Poincare stroboscopic points we can see a general view of a strange chaotic attractor. On the 
other hand Fig. [2b shows system behaviour state for a smaller excitation amplitude (;U = 0.25). 
In this case the initial trajectories follow the strange attractor complex geometry but after long 
enough time {t = 1200) collapse to a steady state motion. Here the corresponding Lyapunov 
exponent Ai = —0.023 is negative which means a convergent motion. 

For better clarity we also plotted the corresponding time histories (Figs. and b) using the 
same colour as in Figs. EK and b. One can easy note that for /i = 0.27 (Figs. [5^ andEK) we have 
chaotic motion of the system while for /i = 0.25 (Figs. Eb andEb) after some time of chaotic 
transient motion with a large amplitude of oscillations and jumping between potential wells we 
got steady state regular vibrations of much smaller amplitude located in a single potential well. 



3 Summary and Conclusions 



We have examined criteria for transition to chaotic vibrations in the Duffing system with a 
damping term dpt(f ) = described by a term fractional exponent p. We were interested 

especially in p < 1.0 {p = 0.1) representing an important case which can be associated with a dry 
friction phenomenon |22j . The critical value of excitation amplitude ^ above which the system 
vibrate chaotically has been estimated, in the first step, by means of the Melnikov method and 
later confirmed by calculating corresponding Lyapunov exponent. 

The Melnikov method, is sensitive to a global homoclinic bifurcation and gives a necessary 
condition for excitation amplitude /U = /ici system in its transition to chaos [HI E]. On the 
other hand the largest Lyapunov exponent ,23.^, measuring the local exponential divergences 
of particular phase portrait trajectories gives a sufficient condition ^ = for this transition 
which has obviously a higher value of the excitation amplitude /i = > /^ci • 

Above the Melnikov transition predictions {ji > fici) we have got transient chaotic vibrations 
[3 IHl El as we expected drifting to a regular steady state away the fractal attraction regions 
separation boundary. This is typical behaviour of the system which undergoes global homoclinic 
bifurcation. In the region of resonance the amplitude of vibration depends on the damping 
exponent in a nontrivial way JS|- In particalar amplitude response rised for small fractional 
exponents (p < 1). Namely, such increase can be identified as an analogues of increasing system 
response on stronger external forcing. Thus one may expect that a chaotic threshold would be 
reached easier in that case. In fact, comparing Figs. and b, we observe larger Lyapunov 
exponents for a fractional p (p = 0.1 in Fig. however from our Melnikov theory analysis we 
draw a different conclusion. The cases for p < 1 (Figs. Ot'-b) seem to be more stable against 
transition to chaos then the p = 1 case. This interesting contradiction should be clearyfied by 
further simultaneous studies of local bifurcations and resonances in the above system. 

For small p and in the limit of a small velocity we have found a region of system parameters 
leading to an interesting dry friction effect. Note, in our approach this limit has been obtained 
naturally on the same footing as other realizations. This enabled us to compare results of 
various damping terms (Figs, andEt"). However to get reliable results in that region one 
must improve a simple numerical algorithm. We have not studied this effect systematically in 
this paper leaving it for a future publication. 

The present approach can be used to generalize models of magnetorheological dampers in 
novel studies of their inflence on vehicle dynamics [T^. To reduce harmful vibrations one can 
consider application of dampers composed of fractional terms. 
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Figure 1: External potential V{x) for 6 = —1 and 7 = 1 (Fig. and velocity (v = x) 

dependence of a damping term dpt(x)/a for different p normalized by a (Fig. ^p). 




Figure 3: Critical value of excitation fic versus to (Fig. ^) and exponent p (Fig. ISb) for 5 = — 1, 
7 = 1. Points in Fig. shows numerical results of transitions to chaos obtained by Lyapunov 
exponents. 
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Figure 4: Lyapunov exponents for (5 = —1, 7 = 1 and p = 0.1 (Fig. p = 1.0 (Fig. EJ)). 





Figure 5: Phase portraits and Poincare maps for 6 = —1, 7 = 1 and /x = 0.27 (Fig. ^^), 
/X = 0.25 (Fig. Eb). 



